Fully self-consistent LDA+DMFT Calculation: A plane wave and projector 
augmented wave implementation and application to Cerium, Ce203 and PU2O3. 
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'CEA, DAM, DIP, F 91297 Arpajon, France 

The combination of density functional theory in the local density approximation (LDA) and dy- 
namical mean field theory (DMFT) has been successful to describe localized or delocalized correlated 
electrons in condensed matter physics. However, the accurate calculations of structural properties in 
this framework are scarce and drastic simplifications are sometimes made, such as the atomic sphere 
calculation (ASA) or the lack of self-consistency over electronic density. In this paper, we present 
an implementation of the full self-consistency over electronic density in a projected local orbital 
LDA-I-DMFT framework on the basis of a plane wave-projector augmented wave DFT code. This 
allows for an accurate calculation of total energy. We show the application to Cerium, Cerium oxide 
Ce2 03 and Plutonium oxide PU2O3. We emphasize the improvement due to both the self-consistency 
and the precise description of the density (i.e. without the ASA). In order to have a correct and 
physical calculation of the energy terms, we find that the calculation of the self-consistent density 
is mandatory. We compare LDA-I-DMFT and LDA-f U solutions, and underline the qualitative dif- 
ferences of their converged densities. In particular and as a consequence, LDA-I-DMFT implies, 
on top of a better physical description of correlated metals and insulators, a reduced occurence of 
unphysical metastable solutions in correlated insulators in comparison to LDA-I-U. 

PACS numbers: 71.15.-m,71.30.-f h,71.15.Ap,71.15.Mb 



I. INTRODUCTION 

Electrons in localized orbitals exhibit strong inter- 
action effects. Their ab-initio description, in con- 
densed matter physics, have made progress thanks to 
the developpement of the LDA+\]^ and LDA-^DMFT 
methodsS^^S. This last method has in particular been very 
successful to describe both the itinerant and localized be- 
havior of strongly interacting electrons. 

Whereas LDA-I-U is implemented routinely in DFT 
electronic structure codes to compute all kinds of prop- 
erties, frameworks to implement LDA-I-DMFT for an ar- 
bitrary basis (e.g plane waves) have emerged only re- 
cently. Two different schemes both using Wannier func- 
tions were used: the construction- of Maximally Local- 
ized Wannier Function (MLWF)^ or NMTO^o and the 
calculation of Projected Local Orbitals (PLO)ii"— . The 
PLO framework is light and does not require the cal- 
culation of MLWF: In particular, MLWF for f-electrons 
systems might be technically difficult to compute. 

Besides, the need for a precise method to compute 
total energy requires the self-consistency over the elec- 
tronic density to be performedi^"— . Calculations of 
total energy using the self-consistency over density are 
scarcei^^ii^. 

In this paper, we present an implementation of 
LDA-I-DMFT with total energy in the PLO scheme, using 
the self-consistency over density in a plane wave (PW) 
based PAW code. This implementation is done in the 
code ABINIT—"^ and is an extension of a previously 
published non self-consistent (NSCF) projected local or- 
bital scheme in PAW—. Since then, a similar NSCF im- 
plementation from the output of the PAW code VASP^ 
has been described^^. The Projector Augmented Wave^l 
framework in combination with e.g a plane wave basis is 



a widely used tool^*"^° to carry out accurate electronic 
structure calculations. The framework of our implemen- 
tation is general and could also be used in e.g real-space 
based PAW code^i^. Moreover, our implementation put 
the computation of energy in LDA-I-DMFT on the same 
footing as LDA : a wide range of systems are available 
for description. 

The plan of the paper is the following. In Sec|lll the 
PLO framework for DMFT is briefly described, expres- 
sions for total energy are given and the way the self- 
consistency is performed is detailed. We emphasize in 
particular that having two expressions for the total en- 
ergy is a check to avoid any errors in the implementation. 
Technicalities related to the PAW formalism are given in 
Appendix |^ for self-consistency and Appendix [B] for the 
PLO scheme. 

In Sec lIIIl we study three correlated systems with 
/ orbitals. For these systems, our PLO scheme for 
LDA-fDMFT is particularly adapted: using MLWF 
would be possible, but the determination of MLWF is 
presently not a straightforward task for / electrons sys- 
tems. We first investigate cerium and cerium oxide 
Ce203. A previous implementation of self-consistency in 
LDA-f DMFT using an ASA basis was applied to these 
compounds. This is thus the opportunity to put forward 
the improvement brought by our approach on these sys- 
tems: because of the PW-PAW basis, the physical accu- 
racy is better, and the method is more flexible. Besides, 
we illustrate, for the sake of clarity, the difference be- 
tween self-consistent densities obtained in LDA+U and 
LDA+DMFT for CcaOa. We then present a study of 
PU2O3. We compare this compound to cerium oxyde 
and our results to previous LDA+U calculations. For 
these systems, we emphasize that LDA+DMFT should, 
at least in the case of cerium and cerium oxyde, de- 
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crease the number of spurious metastables states as found 
in LDA+U'^'^"'^^. Especially for cerium, on the basis of 
our calculations on 7 cerium, the calculation of complex 
phases of cerium, such as /3 cerium, should be easier in 
comparison to LDA+U'^^ . 

In Appendix [B1 we show that our implementation can 
be precisely compared with an DFT-LDA+U implemen- 
tation, when the self-energy in DMFT is chosen to be the 
static mean field self-energy. We illustrate this on the 
case of cerium. The agreement between the two calcula- 
tions shows the consistency of the two implementations. 
In Appendix [Cl we discuss the impact of the windows of 
energy used to compute Wannier orbitals on the spectral 
function of cerium: this underline the need for a con- 
sistent determination of U with the choice of the energy 
window. 



II. THEORETICAL FRAMEWORK 

A. LDA+DMFT formalism: Projected local 
orbitals 

We summarize here the main equations from the Pro- 
jected local orbitals framework^^iir— . In this framework, 
a subspace W of the total Hilbert space is used as a basis 
for the DMFT calculation. This subspace is spanned by 
a given number of DFT Kohn Sham (KS) Bloch wave 
functions. The local orbital on which the coulomb inter- 
action applies are defined on the following way; The KS 
Bloch wave function ^Pkiy are projected on local orbitals 
Xkm (*5-§- atomic orbitals). Let's define P,-^^(k). 

p^.(M) ^ (xLl*k.) (ii.i) 

where m = 1 , . . . , M is an orbital index within the cor- 
related subset, R denotes the correlated atom within the 
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B. LDA+DMFT formalism: Self-consistency over 
density 



The scheme of the self-consistency over electronic den- 
sity within LDA-I-DMFT is shown on Fig. [1] and briefly 
described in the caption. In more details, from the 
DMFT loop, the non-diagonal occupation in the Bloch 
band index can be extracted from III. 5] and used to com- 



primitive unit cell, k is a k-point in the Brillouin Zone, 
and v is the band index. In Appendix IB 21 we discuss 
the calculation of the projection in PAW. This projec- 
tion defines Wannier functions. These Wannier functions 
are built upon the local orbitals and are an orthonormal- 
ized linear combination of a limited number of KS wave 
functions. Let's define: 

IXkJ^ E<*k.|xDl*k.) (II.2) 

The orthonormalization of |x^„j) gives the Wannier func- 
tions lui^). We call Fj^^(k) the projections accordingly 
orthonormalized . 

The flexibility in the choice of W, and thus in the 
choice of Wannier functions enable us an easy compari- 
son to calculations using Maximally Localized Wannier 
Functions^iii^ or atomic orbitals (which are used in the 
PAW LDA-HU implementations^iS). 

From the orthonormalized projections, the local quan- 
tities can be linked to the quantity deflned on the lat- 
tice. For the Green's function and the self-energy, we 
thus have: 

G W ' i^^n) = E (k) Gil (k, tUJn) P^*m' (k) , 

k,(i/iy')e>v 

(II.3) 

AS;^|„(k,za.„) = E E ^.™(k)AS3(*^n)P,^v'(k), 

R. mm' 

(II.4) 

where 



(11.5) 
(11.6) 

I 

pute the density using: 

n(r) = (r|n|r) = E (r|*k^)(*k..|fi|*k,v')(*k^'|r) 
= E (r|*k.)/...'.k(*k.|r) (II.7) 

k,i/,i^' 

Then, the DFT Hamiltonian is built and diagonalized: 
New KS eigenvalues and eigenfunctions are extracted. A 
special care is taken to obtain the new KS bands with a 



[ [{iuJn + Ai - £ku)Su^' - ASbi(k, iwn)] ^ \ 
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FIG. 1: Fully Self-consistent LDA+DMFT scheme adapted 
to the Projected Local Orbital scheme used in our implemen- 
tation. After each solution of the DMFT loop, non diago- 
nal occupations in the Kohn Sham Bloch basis are used to 
compute the total electronic density (in practice, its PW and 
PAW components). This electronic density is used to build 
the Hamiltonian which is diagonalized to give new KS wave 
functions, and eigenvalues. 



good accuracy for each new electronic density, especially 
for unoccupied KS states. Then, proicctions III. II are re- 
computed to build Wannier functions for the next DMFT 
loop. KS eigenvalues are also used to compute the Green 
function using ilLSl 

Pecuharities of the PAW formalism are described in 
Appendix |Al Following the discussion of Ref.— , we thus 
confirm that our implementation allows the evaluation of 
the self-consistent charge. 



{Hu) and -Bdc depend only on local quantities and are 
computed from the resolution of the impurity problem. 
With the HI solver, {Hu) is computed with the Migdal 
formula with a scheme similar to the one detailled in 
Refii.^. i?LDA['T'DMFT(r)] is computed from the knowl- 
edge of the density computed in Eq. III. 71 (see Appendix 
|D|) . The full localized limit (FLL) version of the double 
counting is used here as in RefJ^. 

The kinetic energy can be computed from the one par- 
ticule non diagonal density matrix ni(r,r') — (r|n|r') as: 



1 d^ni{r' — r, r) 



dr 



= - dr /^.<^'.k*,..k(r) — *^'.k(r) 

= E f^^^t^.^'.^ (11-9) 



The band energies can be computed as: 



C. Calculation of Internal Energy in LDA+DMFT 

The LDA+DMFT formalism can be derived from a 
functional of both the local density and the local Green's 
function^. The internal energy can be derived and 
writes^: 

-£;Ha[n(r)]-|-Bxc[n(r)]-/ v^,n{r)dr 

, " s 

^'lDA+DMFT = £^LDA['^DMFT(r)] + Tr[i/LDAGida] 

+Tr[HLDAGdmft] + {Hu) - Euc 

(II.8) 



Tr[GLDA-ffLDA] — Y^ f^^^[nuMFT]£li^^[nuMFT] 
Tr[GDMFT^^LDA] = /,'?k'^'^ [»DMFT]e|^,k^ [»DMFt] 



Thus we have the two following expressions for the 
total energy in LDA+DMFT: 



^LDA+DMFT 



E /.T^^f^.k + i^LDA BCHV)] +{Hu) ~ Euc 



iy,k 



-£;Ha[n(r)]+Bxc[n(r)]-/ i>xc"(r)dr 



(11.10) 



^^£da+dmft = E /°^r*-^-'^k + i^xc+Ha[ri(r)] + / dr«ext(r)n(r) + (ffc;) - Euc (11.11) 



The first expression is the so called double counting ex- 
pression and the second one is the direct one. As in the 
usual DFT calculations in ABINIT^i^, the exact nu- 
merical equality of these two expressions is used in our 
calculations, both as a test of correctness of the numerical 



scheme and of the convergence. 



This implementation has been done in the code 
ABINITi2ai2^ using the PW-PAW implementation^. 



4 



III. APPLICATIONS 

We show here the apphcation of the self-consistent 
LDA+DMFT implementation to Cerium, Ce203 and 
PU2O3. The first two systems have been already 
studied with a fully self-consistent LDA(ASA)-HDMFT 
technique^*. This implementation was based on an 
LMTO-ASAii framework. We will thus be able both to 
check if physical effects due to self-consistency are recov- 
ered in our implementation and what are the improve- 
ment brought by the use of the PW-PAW method. In- 
deed, the PW-PAW implementation is more flexible and 
more precise: the inclusion of semi-core states is easier 
and structural properties are expected to be more pre- 
cisely computed. 

We have performed both LDA-I-U calculations and 
LDA-fDMFT calculations. Two solver were used for the 
LDA-I-DMFT calculations. The first one - the static 
Hartree Fock solver - , is used in order to check that 
LDA-I-U results are recovered with this approach (see 
Appendix IB 31) and thus to check a part of the imple- 
mentation of the self-consistency over electronic density. 
The second solver we used is the Hubbard I solver, as in 
Ref^® (see appendix [E| . For simplicity, we restrict our 
LDA-f DMFT calculations to the paramagnetic phase of 
these compounds. 



A. Cerium 

Cerium metal exhibits an isostructural phase transi- 
tion between the large volume (7) phase and a small 
volume (a) phase. In the 7 phase, electrons are local- 
ized, whereas upon pressure, electrons and more and 
more delocalized^. Several studies using a variety 
of methods have been used to understand the elec- 
tronic structure of the 7 phase or the mechanism of the 
transition,^^ii^i^""^ In this LDA+DMFT study, we are 
using a simplified solver for the impurity model. We thus 
study only the 7 phase in order to compare our results 
to the calculations of Pourovski et aA^. We focus here 
on the spectral function and on structural properties. 
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FIG. 2: Spectral function of 7-cerium, computed in 
LDA+DMFT (HI) with and without Self-consistency on the 
charge density 
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FIG. 3: Energy versus volume curve for 7 Cerium computed 
in the LDA+DMFT framework with the Hubbard I solver. 
As discussed in the text, this energy variation comes mainly 
from -Elda[wlda+dmft]. 



Appendix [C|) . Calculation of the spectral function is car- 
ried out at the experimental volume (34.8 A'^). 



2. Spectral functions and structural parameters: Results 
and discussion 



1. Computational scheme and details 

In the PAW datasets, 5s,5p,6s,5d and 4f states are 
taken into the valence. Atomic data for PAW are taken 
from Ref.'^^. The cutoff energy for the plane wave expan- 
sion is 12 Ha. We found that having 28 k-points in the 
Irreducible Brillouin Zone is sufficient to have a preci- 
sion of 0.02 ua on lattice parameter and 0.2GPa on bulk 
modulus. All calculations are carried at 273K. For the 
DMFT loop, 20 KS functions are used to define corre- 
lated orbitals. They constitute the basis for the DMFT 
loop and such as choice should be appropriate to do a 
comparison with ASA calculations (see the discussion in 



Spectral functions are gathered on Fig. [51 and struc- 
tural properties are gathered in Tab. [T] 

First, the effect of self-consistency over density shifts 
the position of both the lower and upper Hubbard band 
towards higher energies by 0.7 eV. In comparison the 
shift was 0.2 eV in the ASA calculation^®. Concerning 
structural properties, the inclusion of self-consistency in- 
duces a change of 2% on the lattice parameters, whereas 
in ASA, the change is only 0.5 %. So the magnitude of 
the effect of self-consistency is different in LMTO-ASA 
and in PW-PAW. Effect of self-consistency over spectra 
and structural properties are thus qualitatively similar to 
the study of Pourovskii et al but the magnitude of the 
effect is larger in PW-PAW. 
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a (a.u.) 


Bo (GPa) 




9.76 


19/21 


PAW/LDA+U 


9.58 


32 


PAW/LDA+DMFT NSCF (H-I) 


9.41 


38 


PAW/LDA+DMFT SCF (H-I) 


9.58 


36 


ASA/LDA+Ui^ 


9.44 


49 


ASA/LDA+DMFT NSCF (H-I)i^ 


9.28 


50 


ASA/LDA+DMFT SCF (H-I)i^ 


9.31 


48 



TABLE 1; Lattice parameter a and Bulk modulus Bo of 
7 Cerium according to experimental data and calculations 
within LDA+U and LDA+DMFT, in PAW and in ASA^^ 
For a appropriate comparison, PAW/LDA+U are carried out 
with the expression of the density matrix shown in appendix 

El 

We then compare ASA and PAW results: the differ- 
ence of absolute values of the lattice parameter found in 
LDA-f DMFT in the two basis is of the same magnitude 
as the difference found in LDA-|-UiS,. In our calculations 
with the parameters used, the agreement of both spectra 
and structural properties with experimen t^^i^^ is slightly 
better compared to ASA calculations^^. 



3. Origin of the variation of the energy versus volume 

We discuss now the behavior of the energy versus vol- 
ume curve in LDA-I-DMFT. We can split the energy in 
two contributions: the first one include the interaction U, 
namely {Hu) — E^c and the second include all the other 
terms and is i?LDA ['T'LDA-i-dmft] • If we separate these two 
contributions as a function of volume, one thus see that 
the first term is negligeable for Cerium, in LDA-I-DMFT 
(Hubbard I). It comes from the fact that the number 
of f-electrons found in the solver is approximatively one 
in this approximation, for each volume. In other words, 
the variation of energy as a function of volume is com- 
pletely described by i?LDA [?^lda+dmft] • It means that, 
for each volume, LDA-I-DMFT converges to a new den- 
sity, riLDA-i-DMFT and that the LDA energy alone for this 
set of densities show the behavior plotted on FiglH 

In LDA+U, the energy decomposition show a slightly 
different behavior: the term (Hu) — i?DC is not com- 
pletely negligeable. In fact, computing the energy vari- 
ation with only -Elda [?^lda+u] leads to a lattice pa- 
rameter of 9.35 ua (about 2% less than the LDA+U 
value). It comes directly from the fact that the num- 
ber of correlated electrons in LDA+U changes as the 
volume changes and modify (Hu) — Edc- However, the 
main effect is the same in LDA+U and LDA+DMFT: the 
change of the electronic density induced by correlation, 
shows a variation as a function of volume, which shifts 
the equilibrium volume mainly because of the variation 
of -EldaI^lda-i-dmft]. 

Let's now analyse the variation of the converged den- 



sities as a function of volume with respect to the en- 
ergy. In fact, we observe that the variation of i?LDA['T'] 
for a given density, can be linked to the number of f- 
electrons: for a variation of volume of 12% around the 
equilibrium volume, the number of electrons increases of 
0.007 in LDA, and decreases of 0.002 in LDA+U and 
of 0.03 in LDA+DMFT . Clearly, the Hartree interac- 
tion energy should thus decrease as a function of volume 
in the case of LDA+U/LDA+DMFT. In all these three 
methods, however the impact of these variations are not 
the same on the interaction energies. If might come from 
the different density matrices which are the results of 
these three calculations: In LDA+DMFT, electrons are 
more homogeneous distributed because of possible fluc- 
tuations between different localized atomic states. On 
the contrary, in LDA+U, the electron is localized in a 
single orbital, whereas in LDA, electrons are distributed 
among orbitals according to crystal field. In other words, 
Hartree interactions between electrons belonging to dif- 
ferent orbitals are not the same and should be lower than 
the Hartree self-interaction of one electron in one orbital. 
As a consequence, the variation of the number of electron 
as a function a volume is a hint to understand the varia- 
tion of energy, but a complete understanding requires the 
knowledge of the density matrix, or of the total electronic 
density. 

In the case of non-self consistent LDA+DMFT calcu- 
lations, the expression lll.lOl is used for the total energy^. 
The main effect comes from the band energy. The dif- 
ference between the DMFT band energy and the LDA 
band energy decreases when the volume increases as pre- 
viously observed^^. However, the intensity of this vari- 
ation is weaker than the correct variation at the self- 
consistent density. Moreover in the SCF calculation, this 
large variation is partially compensated by LDA double 
counting terms. It show that to have a physical correct 
behavior of the different terms of the energy, achieving 
convergency to the LDA+DMFT density is mandatory. 
As emphasized before, this is the very change of the elec- 
tronic density which make the difference in i?LDA- Such 
change cannot be reproduced by a non-self-consistent cal- 
culation. 



4- Search for the ground state of the system, metastable 
states 

As underlined in the previous paragraph, the or- 
bital anisotropy obtained in the ground state with 
LDA+DMFT is very small and is thus different from 
the ground state obtained in LDA+U which has a large 
orbital anisotropy^^^. We will illustrate the same re- 
mark in the case of cerium oxyde in the next subsection 
with plots of the electronic densities. We thus expect the 
convergency in the case of complex system with large 
unit cell (such as /3 cerium) to be much easier with the 
LDA+DMFT method than with the LDA+U method^. 



6 



B. CeaOa 

Ce2 03 is a Mott antiferromagnetic insulator with an 
optical band gap of 2.4 eV— and a Neel temperature of 
9K. A ionic counting of electron leads to one / electron 
in the valence band. The qualitative picture of an insu- 
lator has been obtained by first principles calculation us- 
ing the LDA+\&r^, Hybrid functional^, LDA-hDMFT 
frameworki^, and GW-|-LDA-|-U^. In most all these cal- 
culations, one electron is localized and the ground state 
is an antiferromagnet. LDA-I-DMFT is able to describe 
a paramagnetic insulating solutionis. We present here 
the application of our LDA-I-DMFT scheme to this com- 
pound using the PW-PAW framework. 



1. Computational scheme and details 
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FIG. 4: Spectral function of Ce2 03, in LDA-^DMFT (Hub- 
bard I) with and without self-consistency. 



In this calculation, we do not relax the internal param- 
eters: the volume is the only variable of the system. For 
Cerium, we use the same atomic data as above. For oxy- 
gen, 2s and 2p electrons are in the valence. PAW match- 
ing radius for oxygen is 1.52 ua. The cuttoff energy for 
the plane wave expansion is 30 Ha. We use 32 k-points in 
the Irreducible Brillouin Zone. Structural properties are 
converged to less than 0.01 ua for the lattice parameter 
and less that 0.1 Mbar for the bulk modulus. 

Calculations are carried out at a temperature of 273 
K. 52 Kohn Sham bands are used to define correlated 
orbitals and thus as a basis for DMFT calculations. As 
for cerium, this number of bands should be appropriate 
to do a comparison with previous calculations with the 
ASA formalism using a basis containing the same number 
of bands. The value of U is taken to be 6.0 eV and J is 
eV. 



a (A) ^0 (Mbar) 





3.89 


1.11 


PAW/LDA-hU ( AFM)S 


3.87 


1.3 


PAW/LDA-hU(AFM) 


3.85 


1.5 


PAW/LDA-(-DMFT (H-1) NSCF 


3.76 


1.7 


PAW/LDA-hDMFT (H-I) SCF 


3.83 


1.6 


ASA/LDA -hUi^ 


3.84 


1.5 


ASA/LDA+DMFT(H-I) NSCPi^ 


3.79 


1.6 


ASA/LDA+DMFT(H-I) SCF^^ 


3.81 


1.6 



TABLE 2: Lattice parameter a and Bulk modulus Bo of 
Ce2 03 according to experimental data and calculations with 
various frameworks. In our calculation and in Ref.— , the c/a 
ratio is fixed to the experimental value (1.56). PAW/LDA-I-U 
are carried out with the expression of the density matrix cor- 
responding to Eg IB. 31 and with U=6 eV. Calculations are 
done at a temperature of 273K. Note that entropy is ne- 
glected in these calculations. Calculations from Da Silvai^ 
and Pourovskii et aJi^ use respectively U=5.3 eV and 5.4 eV. 



2. Results and discussion 

Spectral functions are given on Fig. 2] and structural 
properties are gathered on Tab. 

Concerning first LDA-I-U results for the structural pa- 
rameters, our calculations for the ground state, with anti- 
ferromagnetism (AFM) are in good agreement with cal- 
culations of Ref^-^. The slight disagreement concerning 
the bulk modulus might come from the absence of relax- 
ation of internal parameters in our calculations. 

Results from LDA-I-DMFT calculations reveal that ef- 
fect of self-consistency over spectra and structural prop- 
erties are qualitatively similar to the study of Pourovskii 
et al: The volume increases and the gap slightly decreases 
in the self-consistent calculation with respect to the non- 
self-consistent one. The quantitative results are however 
differents from Refii^: the ASA approximation induces 
an slight error of 1% in the determination of the lat- 
tice parameter with respect to full potential codes. We 
correct this error in the PW-PAW implementation of 
LDA+DMFT calculations. 

Spectral function computed with LDA-I-DMFT with 
the PW-PAW basis show a good agreement with the op- 
tical gap of 2.4 eV— as encountered also in LDA-|-U^. 



3. Comparison of LDA, LDA + U and LDA+DMFT 
densities 

The analysis of the variation of different energy terms 
as a function of volume leads to the same conclusion as for 
Cerium: The variation of energy is mainly due to -Blda 
and is correlated to the variation of the number of / 
electron. As the volume increases, this number increases 
in LDA, decreases in LDA-I-DMFT and is nearly constant 
in LDA-f U. 

Also, as observed in the case of cerium, the orbital 
anisotropics are different in the LDA-I-U results and the 
DMFT calculation. In particular, one f-orbital is filled 






FIG. 5: Difference between electronic densities computed in 
the LDA+U approximation and in the LDA approximation 
viewed along different direction of the cell of Ce2 03. The 
density of the filled orbital f^~ is clearly visible. Blue (resp. 
green-red) area corresponds to positive (resp negative) value 
of the difference. 



with the LDA+U method whereas in DMFT, electrons 
are allowed to fluctuate between different orbitals, hence 
the diagonal terms in the occupation matrix of the local 
orbital basis are similar and not far from 0.07. 

To illustrate the differences between LDA+DMFT and 
LDA+U electronic densities, we have plotted^ on Figl5] 
(resp. Fig. [H]), the difference between LDA+U (resp. 
LDA+DMFT) and LDA electronic densities, along dif- 
ferent crystallographic directions. In agreement with the 
density matrix, LDA+U shows a localization of one elec- 
tron in only one f-orbital and in LDA+DMFT the f- 
electrons are spread among f-orbitals, however, in a dis- 
tinct way, with respect to LDA. This difference of be- 
havior is not visible on the total spectral function which 
contains the sums of the contributions of all orbitals. A 
plot of the density of states resolved in the projection of 
the orbital moment would show the same kind of infor- 
mation. 



4- Thermodynamical study 



FIG. 6: Difference between electronic densities computed in 
the LDA+DMFT approximation and in the LDA approxi- 
mation viewed along different direction of the cell of Ce2 03. 
Contrary to the case of LDA+U, one encounter in this case 
only a slight redistribution of electrons inside f-orbitals. Blue 
(resp. green-red) area corresponds to positive (resp negative) 
value of the difference. 



We have thus minimized the structure of the molecule 
O2. The distance is found to be 1.22 A'^ and the cohesive 
energy is 7.6 eV in agreement with previous results (see 
e.g. Ref..^). 

Aril computed in LDA+DMFT with the Hubbard I 
solver gives thus a value of -19 eV ± 0.1 eV, within 1. eV 
of the experimental results (-18.63 eV)^ and in the range 
of existing value computed in LDA, and LDA+U^"— . 
Upon variation of U between 5 and 7 eV, the computed 
energy show a change of less than 0.15 eV. 



C. PU2O3 

PU2O3 is an paramagnetic metal insulator above 10 
K— . Its conductivity gap is 1.8 eV and in a ionic pic- 
ture, it should contain 5 / electrons. In this section, we 
describe this compound with GGA(PBE)+DMFT (Hub- 
bard I) and compute the spectral properties of this sys- 
tem and we discuss the improvement with respect to 
GGA(PBE)+U36.70. 



From the energies obtained for cerium and cerium 
oxyde, one can compute the variation of internal energy 
Aril of the following reaction: 



2Ce+ -O2 ^ CesOa 



(III.12) 



1. Computational scheme and details 

We use the same parameters as for cerium oxyde. For 
the PAW basis, semicore states are included in the atomic 
dataset'^^ and for the PW basis, the cutoff energy for the 
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plane wave is 30 Ha. The value of U is taken to be 4 eV. 
As for cerium oxyde, 52 Kohn Sham bands are used to 
compute Wannier / orbitals. 

2. Results and discussion 

Spectral functions for the SCF and NSCF calculations 
are shown on Fig \7\ The GGA+DMFT calculations de- 
scribes PU2O3 as a Mott Hubbard insulator with a f-f gap 
of 1.8 eV. This is similar to what is obtained with the 
GGA+U method^ and near the experimental conduc- 
tivity gap of 1.8 eV— - which contains, however two par- 
ticles excitations. As for cerium oxyde, self-consistency 
shifts the f peak towards higher energy with respect to 
the 0-p band and to the conduction band. 

The density matrix obtained in the paramagnetic so- 
lution is different from the LDA-I-U results. However, 
contrary to the case of Ce203, there is still an orbital 
anisotropy coming mainly from crystal field : crystal 
field splitting is indeed larger in plutonium oxyde, in 
agreement with the fact that Cerium / orbitals are ex- 
pected to be more localized than Pu / orbitals. From the 
energy levels computed in Hubbard I (expression lE.ip , 
we can estimate the crystal field splitting to 0.05 eV in 
Ce203 and 0.5 eV in PU2O3. It thus explains why or- 
bital anisotropy is observed in PU2O3 and not in Ce203: 
It comes from the direct physical effect that lower levels 
have more statistical weights and in PU2O3, the difference 
in energy between levels is sufficient. Indeed and logi- 
cally, it is the same orbitals that are filled in LDA-|-U2& 
that are partially filled in LDA-I-DMFT. Nevertheless, in 
LDA-I-U, one can encounter spurious orbital anisotropy 
and associated metastables states that are due to a break- 
ing of symmetry, or to the filling of higher states in 
energy and the subsequent stabilization of their energy 
levels. Such breaking of symmetry does not happen in 
LDA-I-DMFTS, because this method takes into account 
all possible slater determinant in the Green's function: 
metastables states encountered in LDA-f U which have a 
nearly degenerate energy are treated on the same foot- 
ing in the DMFT Green's function. However we cannot 
exclude metastables states totally in the case of a large 
crystal field or spin-orbit coupling. However they should 
be highly metastable in these cases. 

IV. CONCLUSION 

In this paper, we first present a self-consistent imple- 
mentation of LDA-I-DMFT using the projected local or- 
bital formulation. This implementation is done inside 
a plane wave-projector augmented wave code. We de- 
tail the PAW specific implementation and two general 
expressions of the total energy. 

Then, we apply this method to three strongly corre- 
lated systems with / orbitals. For these systems, the 
PLC scheme is particularly adapted, because building 
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FIG. 7: Spectral function of PU2O3, in GGA-hDMFT (Hub- 
bard I) with and without self-consistency. 

MLWF for / systems is not a straightforward task. First, 
we study cerium and Ce203. We highlight on these sys- 
tems the improvement brought by the plane wave and 
projector augmented wave formulation with respect to 
ASA: The accuracy of the basis is better. Therefore, 
structural parameters (lattice parameters and elastic 
properties) are improved. More generally, with this im- 
plementation, the calculation of energy in LDA+DMFT 
is on the same footing as the same calculation in LDA: 
semicore states are easily included, energy of different 
structures can be easily compared. We underline the im- 
pact of the definition of local orbitals - and thus of the 
energy window - on spectral function and structural pa- 
rameters in the case of cerium. 

Finally, we study plutonium oxyde. We compare this 
compound to cerium oxyde and underline the reduced 
occurence of spurious metastable states in LDA-I-DMFT. 

Therefore, LDA+DMFT, besides giving a better phys- 
ical description of correlated systems, could thus solve at 
least partially the problem of multiple minima which pre- 
vents to find the ground state of large systems with the 
LDA+U methodi^"— . The computational cost could be 
much reduced by using a simplified solver, such as Hub- 
bard I, provided the system studied is a Mott insulator 
and not a correlated metal. To be predictive, however, 
a quantum Monte Carlo solver should however be used. 
The interface between quantum Monte Carlo solvers and 
our PLC implementation of DMFT is in progress. 

The general LDA+DMFT framework with the accu- 
rate PW-PAW basis open the way to the calculations of 
others properties, such as forces, phonons^. Our im- 
plementation has been done inside the open source code 
ABINIT^-^. 
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implementations are quite different. In this appendix, 
we briefly compare the definition of local orbitals in these 
two methods, we show how we do it in the PAW method, 
and we give a numerical careful comparison with our im- 
plementations. 



Appendix A: Self-consistency over charge density: 
peculiarities of the PAW formalism 



From InTfl we have (we neglect spin in this section): 
n(r) - ^,^''.k(*.,k|r)(r|*,.,k> (A.l) 

Then, using relations (11) of Reff^ for the operator 
r)(r|, and Ea lA.ll we have: 

nir) = J2 ^,^'Mc*:,k(r)$.',k(r) + 

J2 /-,-Mc(*.,k|p*)(((p..|r)(r|(/j,) - 

ij i/,z^',k 

(^,|r)(r|^,))(p,|$,,,k) (A.2) 

(A.3) 

where pj is the projector number nj angular momenta Ij , 
and its projection nij on atom R. Index for atoms are 
neglected on p, (p and ip. thus: 



n(r) = n'{r) + (r) — (r) 



(A.4) 



with 



n'(r) = Y /...'.k«':,kW*.'.k(r) (A.5) 
ii'(r) = ^p^VP.(r)^,(r) (A.6) 
i^'W = ^p-j(^»(r)^j(r) 



and 



= Y ^,^',k(^^,kbj)(p*|«'^',k) 



(A.7) 



(A., 



1^,1^' ,k 



For the on site terms of the density, we just have to 
compute p'^j, and then replace pij by p[j everywhere in 
the formalism. 



Appendix B: Definition of local orbitals and 
comparison of LDA-|-DMFT and LDA+U 
implementations 

LDA-f U and LDA-f DMFT can be seen as two differ- 
ent approximations for the same hamiltonian. LDA+U 
is a static mean field solution, and LDA-I-DMFT con- 
tains dynamical fluctuations. However, their practical 



1. Definition of local orbital, and density matrix in 
the PLO scheme 



Projections in the PLO scheme are defined by Eq. III. II 
We assume that Xm ^ set of orthonormalized orbitals, 
such as MLWF or atomic orbitals. If we use an infinite 
window for the construction of Ix^n) (Eq |II.2p then we 
have Ixkm) = \Xkm)^ because the basis of Kohn Sham 
functions is complete. In this case, the construction of 
Wannier orbitals from Xm j^st a renormalization with 

Thus the density matrix of correlated electrons writes: 



3 



ix^lx^) 



as it can be rewritten: 



R 

Am,(T 



(B.l) 
(B.2) 



=), (B.3) 



with n the density matrix which is ri = 

Ek,./k%-l*U(*k,.l 

If we assume that x^m ^ orthonormalized 
orbitals then (x^lx^)=l- 



2. Local orbitals for PLO in PAW 

In our PAW implementation of the PLO scheme for 
LDA+DMFT, we follow the lines of the Ref.i^. Equation 
inithus writes as Eq (A.2) of Ref.-i^ in the PAW method: 



nt(k) 



(X^l^^)(?5^l*k.) 



(B.4) 



is the projector n for angular momenta I, and its 
projection m on atom R and Xm is chosen to be the 
atomic eigenfunctions for the given angular momentum. 

The Xm '^ot completely orthonormalized because in 
our implementation, they are used only inside the PAW 
sphere. However, 98% of their density is located inside 
the sphere, so in this discussion, we neglect this slight 
non-orthogonality. 

As we use only the part of the density matrix computed 
inside the sphere with the all-electron wave function (see 
Ref ]^^'^^ ), (x^lx^) represent the integral of the atomic 
wavefunction inside the PAW sphere. We thus have, in 



10 



10 



5 









' 1 ' 


1 ' 

A 




-■- LDA+DMFT(HF) , 


ll 


- 


LDA+U 


h' 
■I 

\.\ 

1 ^"Txy 




..J: 1 



co(eV) 



FIG. 8: Spectral function of 7-cerium, computed in LDA+U 
and LDA+DMFT (HF) 

the PAW framework, with (m index are removed for sim- 
plicity in the equation because integrals are done here 
over radial part of wavefunction) 

Thus, having the same local orbitals in the 
LDA+DMFT and LDA+U frameworks requires to use 
Eg IB.Sl in LDA+U and using a complete KS basis set in 
LDA+DMFT. We show the application of this compari- 
son in the next subsection. 

Another way of computing the pro iection III. 1 1 and oc- 
cupation matrices exist, such as an integrated value in 
PAW sphere of angular-momentum-decomposed charge 
densities^^"^. The slight drawback of such a scheme is 
that the starting atomic orbitals |x^) are not easily de- 
fined in this approach. However, the advantage of this 
approach would be that the spectral weight is by defini- 
tion not lost^ - if the projection is done in all space with 
all terms of Eq. (A.l) of ReL^^. 

As for simplicity, we have decided to do the projec- 
tion only in the sphere, we have used Eqs. IB. 41 and IB. 51 
Differences between the two schemes are smal P^i^^ . 

We underline that atomic data that are used with this 
framework have to be transferable in a large energy range 
because unoccupied excited states are used in the frame- 
work. In particular, the closure relation over KS states do 
not hold if atomic data is not tranferable in the relevant 
energy range. 



ranging from -32 eV to 30 eV (Fermi level is at zero en- 
ergy). To assess the lack of completeness of the Kohn 
Sham basis, one can compute the value of (x^|x^) using 
the closure relation: 

(x^lE[*k..)(*k,.]lx^) (B.6) 

k,i/(T 

In the case of the PAW data that we use, and the bands 
chosen (see caption of Tab. |3]), we find an average value 
of 0.95, whereas the numerical value of (x^lx^) is 0.964. 
The Kohn Sham basis thus describes the atomic wave- 
function with a small error of less that 2%. 

Spectral functions, shown on FigElare nearly identical. 
Structural properties, gathered in Tab. [3]show variations 
smaller than 0.1% on the lattice parameters and 1% on 
the bulk modulus. This is a confirmation of the coherence 
of the two implementations. 





a (a.u.) 


Bo (GPa) 


PAW/LDA+U 


9.58 


32 


PAW/LDA+DMFT(HF) 


9.59 


31 



TABLE 3: Lattice parameter a and Bulk modulus Bq of 7 
Cerium according to calculations using LDA+U with expres- 
sion IB. 51 for the density matrix and LDA+DMFT using 30 
bands ranging from -32 eV to 30 eV around the Fermi level. 



Appendix C: Impact of the energy windows used to 
define Wannier functions on physical properties 

In this section, we study the dependancies of the spec- 
tral function and structural parameters as a function of 
the number of KS states used to define Wannier func- 
tions for the LDA+DMFT calculation. We have carried 
out LDA+DMFT calculations of Cerium for several sets 
of KS states (see Tab. |4|). The choice of the KS basis 
implicity define the local orbital subset as an orthonor- 
malized linear combination of these states. As the locali- 
sation of the local orbital change according to the choice 
of the energy window, a value of the interaction \J should 
be defined for each choice of the set of KS bands. How- 
ever, in these calculations, for testing purpose, we used 
a fixed value of the interaction \J . Spectral functions are 
gathered on Fig. [9] and structural parameters are given 
in Tab. [H 



3. Numerical comparison of the two schemes 

We have carried out this comparison for the volume 
and spectra of 7 Cerium. Parameters for the calcula- 
tion are similar to those of section IlIIAl The calculation 
are carried out at the temperature of 273K. Moreover, 
calculations using LDA+U are done with expression lB.51 
for the density matrix and LDA+DMFT uses 30 bands 



1. Spectral function 

We observe a large variation of the spectral function 
for the different cases studied: In particular lower (resp. 
upper) Hubbard band are shifted from -2.5 eV to -l.eV 
(resp. 3.5 eV to 5 eV) As we use the Hubbard one ap- 
proximation to solve the Anderson model, we can use 
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Eq. IE. II to understand this variation. In this equation, 
the double counting correctioniS is computed with the 
number of electron used in the Hubbard I solver, which 
is nearly one, and does not depend on the window of en- 
ergy. In fact, the shift of the atomic level, comes from 
the definition of the LDA atomic levels (Eq. lE.ip : As 
more and more states Kohn Sham states are included in 
the sum, the mean energy of levels increases tangentially, 
ultimately, until the basis is complete. As the lower Hub- 
bard band visible on the spectrum of Cerium corresponds 
to the removal of one f-electron, it is directly related to 
the atomic levels. Thus, as the energy of levels increases, 
the Hubbard band is shifted. The spectra should "con- 
verge" as the number of KS states increases. However, we 
emphasize that the results of this convergence only cor- 
responds to a given choice of correlated orbitals, namely, 
atomic orbital. Another legitimate choices exists. 

For example, we thus found, that it is especially im- 
portant to use a number of bands equal to 20 per spin, to 
have a good agreement with the study of Pourovskii et 
a^. This could comes from the fact, that in this study, 
the ASA calculation is done with a basis set of also 20 
states per spin (5s, 5p, 6s, 6p, 5d and 4f states) and thus 
is roughly corresponding to the same window of energy. 

A more extensive physical investigation would require 
to compute a different value of U for each window of 
energy (or number of KS states) used to define the cor- 
related Wannier functions: as the windows of energy is 
increased, the orbital are more and more localized. Thus 
the value of U should increase, and the double counting 
term in Eq. IE. II should shift downwards the lower band 
peak. This shift might compensate the upward shift in- 
duced by the first term in Eq. lE.ll 



2. Structural properties 
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FIG. 9: Spectral function of 7-cerium, computed in 
LDA+DMFT (HI) for different choices of the KS basis 



Appendix D: Total energy in LDA+DMFT in the 
PAW formalism 

In this section, we briefiy show that the knowledge of 
fv,vi,vk is sufficient to compute all energy terms in the 
PAW method. We use some notations of Ref.— . Accord- 
ing to the calculation of fv^v'^vk, one may compute pij 
and the pseudo-density with Eq. IA.3I and IA.8I All the 
terms of the energy which depends on the local quanti- 
ties (i?\ E^. E\^, E\^) are computed with p-^, and 

■n} (see Ref42). In the plane wave basis, E^c, £^Ha are 
computed directly from the computed LDA-I-DMFT den- 
sity. The sum over eigenvalues and the kinetic energy are 
directly computed from above ( Eq lll.lOl and Ea lll.lll ) . 
We emphasize that there is no approximation to a tight 
binding model for the computation of the kinetic energy 
or the sum over KS eigenvalues. 



Concerning structural properties, we find that the 
equilibrium lattice parameter increases as the window of 
energy is larger. On the same time, the number of / elec- 
trons is increasing as the number of bands increases. We 
might try to correlate the two evolution. However, we 
stress that a more physical behavior should be expected 
with a more accurate solver like QMC. 



Number of KS bands 


15 


20 


25 


Energy range (Ha) 


1.4 


1.8 


2.0 


Lattice Parameter (ua) 


9.42 


9.58 


9.67 


Number of / electrons 


1.03 


1.116 


1.14 



TABLE 4: Range of energy for the LDA eigenvalues for sev- 
eral sets of KS states. The lowest KS states, corresponding to 
5s orbitals are located -1.25 Ha under the Fermi level. Lat- 
tice parameter and Bulk Modulus computed in LDA-I-DMFT 
with the corresponding Wannier function in the Projected Lo- 
cal Orbital scheme. 



Appendix E: Hubbard I solver and interaction 

To solve the impurity model, we use in this paper the 
Hubbard I method. This method neglects the hybridiza- 
tion between the impurity and the bath in the resolution 
of the Anderson impurity model. We have implemented 
the Hubbard I solver in the same spirit as in Ref^— : we 
impose that the impurity Green's function and the local 
Green's function have the same limit at high frequencies. 
It follows the following definition for the atomic level (see 
alsoM): 

Cn' = E P^:{^KkP^r^*{^) - ^I'c - M- (E.l) 

In this paper, we precise the value of the interaction 
U chosen for each system. We choose J=0 in all calcula- 
tions. 

Lastly, we have found that the Hubbard I approxima- 
tion, used as a solver for the DMFT method can lead to 
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a convergence to spurious non-physical states. In partic- 
ular, a solution of the numerical scheme with a lower 
energy than the true ground states (especially at low 
volume and for small value of U) is found. As empha- 
sized in Ref^'*, due to the limitation of the Hubbard I 
solver, the impurity Green's function is not the same as 
the local Green's function. In particular, without tem- 
perature effects, the number of f-electron is an integer, 
whereas in DMFT, hybridation effects leads to a non in- 
teger number of electrons. In practice however, the num- 
ber of f-electron (in the 7 phase) computed respectively 
from the impurity (in Hubbard I approximation) and lo- 
cal Green's function are respectively 0.99 and 1.10 elec- 
trons. In a spurious solution found in cerium oxyde, the 
numbers are respectively 0.01 and 0.74 electrons. This 
results is clearly inconsistent and should be eliminated 
for physical reasons. The use of a quantum Monte Carlo 
method would suppress this inconsistency, because the 
two Green's functions would be equal. 
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